Active viscoelastic models for cell and tissue mechanics

Living cells are out of equilibrium active materials. Cell-generated forces are transmitted across the cytoskeleton network and to the extracellular environment. These active force interactions shape cellular mechanical behaviour, trigger mechano-sensing, regulate cell adaptation to the microenvironment and can affect disease outcomes. In recent years, the mechanobiology community has witnessed the emergence of many experimental and theoretical approaches to study cells as mechanically active materials. In this review, we highlight recent advancements in incorporating active characteristics of cellular behaviour at different length scales into classic viscoelastic models by either adding an active tension-generating element or adjusting the resting length of an elastic element in the model. Summarizing the two groups of approaches, we will review the formulation and application of these models to understand cellular adaptation mechanisms in response to various types of mechanical stimuli, such as the effect of extracellular matrix properties and external loadings or deformations.


Introduction
Living cells, their surrounding extracellular matrices (ECM) and tissues as a whole exhibit viscoelastic properties, that is, having both an elastic and a viscous behaviour.Mechanical tests are used to characterize the mechanical properties of cells and tissues and help us understand and predict cellular behaviour in healthy or pathological conditions at different time scales [1][2][3][4].For example, experimental studies on cell monolayers devoid of substrate have reported a viscoelastic solid-like behaviour under constant strains [5].In other words, on the time scale of seconds, cell monolayers behave like a viscous fluid and dissipate stress, then reach a plateau in stress on the time scale of minutes, which is a characteristic of elastic solid materials [5].On the other hand, cell aggregates exhibit solid-like behaviour on short time scales, for example, in a few seconds, and fluid-like behaviour at time scales of the order of minutes to hours [4,6,7].In this case, the behaviour of cell aggregates is similar to a viscoelastic fluid material.Additionally, the mechanisms governing the viscous-like behaviour can also differ depending on the time scale of load application [8].For example, on the time scale of tens of seconds to minutes, stress dissipation occurs owing to the turnover of actin filaments and reorganization of the actomyosin network [4,8,9].On the time scale of minutes to hours, stress dissipation mechanisms at the cellular scale, such as oriented cell division and cell rearrangements, can start to influence the response to mechanical loadings [8,[10][11][12].Transitions between fluid-like and solid-like behaviours in living cells can also occur in response to mechanical stimuli.For instance, short-term fluidization has been reported immediately upon strain application, which is usually followed by stiffening [13][14][15][16].
The experimental results are often analysed by mathematical models to capture the important features of the material response, presented in terms of model parameters.These parameters can then be used for classification, comparison and prediction of the mechanical behaviour of cells and tissues subjected to other loading conditions [17].By relating the model parameters to the underlying biological processes on the molecular and cellular scale, the physical meaning of the parameters can sometimes be assigned.Since tissues exhibit time-dependent mechanical behaviours [5,18], a common modelling approach is to consider the cells and tissues as a viscoelastic continuum and to describe their mechanical response from quantitative mechanical interrogations, often stress-strain relationships, in terms of a combination of stiffness and viscosity, or elastic and loss moduli under dynamic loadings.These extracted mechanical properties have long been regarded as disease biomarkers [19,20].For instance, metastatic cancerous cells exhibit lower stiffness than benign cells [21].
Perhaps more importantly, cell mechanics is not only a by-product of the underlying molecular structure but also a means for cells to actively adapt to environmental cues, in service to a preferred cellular function, such as cell migration in wound healing [22] and tissue morphogenesis in development [23].Being viscoelastic in nature, cells can dissipate the imposed stress owing to external strains.However, sometimes this passive response is insufficient to maintain the mechanical integrity of the cell.For instance, cells need to use effective stress relaxation mechanisms such as actin polymerization to prevent tissue fracture.In other instances, under rapid strain applications, cells may need to stiffen or actively pull back to make further deformation difficult, thus preventing further damage to their cytoskeleton.From the adaptation perspective, the evolving elastic and viscous properties owing to the active adjustments of cytoskeleton tension can be considered as a way to facilitate the response and adaptation to external stimuli.Probing active cellular behaviour through the lens of mechanics is particularly intriguing because it offers a window from which the adaptation can be quantitatively examined up close with defined mechanical stimuli.
In the cell cortex as well as at the cell population scale, living matter not only responds to external forces or deformation as any traditional material would but also often exhibits force-generating mechanisms emerging from actin polymerization, adhesion dynamics and actomyosin contractility.This process is referred to as being 'active' in this review.Active matter theory is another continuumbased model used to describe the dynamics of cell cytoskeleton and cell monolayers (reviewed in [9,[24][25][26]).This type of model is based on the theory of liquid crystals and can be employed to describe the mechanics of the actomyosin cortex for time scales longer than the turnover time of actin filaments [27].However, these models rarely capture the way such materials respond to external deformations or stresses; it remains challenging to capture, both experimentally and theoretically, the impact of active processes on the mechanical state of living matter.
Experimentally, a common approach is to subject the cellular materials to tension or compression and then assess their response using stress-strain curves.At the tissue scale, mechanical testing machines can be used to investigate the macroscopic behaviour of tissue samples (figure 1a) [28,29].A similar technique has been developed to study the response of cell monolayers when subjected to compression or tension (figure 1b) [30].In addition, to examine the microscale characteristics of living tissues, atomic force microscopy (AFM) and nanoindentation techniques can be employed, where a probe tip is used to apply forces in the pico-nano Newton range [31,32].
The AFM technique can also be used to study the active and passive behaviour of individual cells [33][34][35] and can be adapted for cell pair studies in single-cell force spectroscopy (SCFS), where a single cell adhered to the cantilever beam serves as the probe tip (figure 1c) [36].Furthermore, various methods are developed to study the response of single cells and cell doublets in tension and compression tests.For instance, in microplate assays, cells are compressed or stretched between a fixed and a deformable plate [37,38], microbead assays use microbeads to apply force via magnetic [39,40] or optical tweezers [41][42][43] (figure 1d), and the micropipette aspiration method involves subjecting cells to negative pressures and measuring their deformation [44][45][46] (figure 1d).This method can also be used for cell pair studies by bringing two cells into contact via micropipettes, that is, dual pipette aspiration (DPA) (figure 1c) [47], or it can be combined with other techniques such as optical tweezers [48].Recently, a new micromanipulator device has been introduced that can directly measure forces in cell pairs under controlled loading conditions, leading to advancements in the precision of the interrogation of cell pair mechanics (figure 1c) [49].Different techniques have varied ranges of resolution and loading rates; thus, one has to consider the application requirements when choosing an experimental technique.For instance, AFM techniques benefit from higher spatial resolution and force sensitivity compared with methods employing micropipette aspiration [50,51].
These experimental techniques are employed to probe cellular response at various temporal and spatial scales.In addition to measuring the global mechanical behaviour of cells and tissues, these techniques can also be employed to examine the mechanical properties of specific components of In cell pair studies, mechanical stimulus is applied through various methods such as a single cell adhered to a tipless atomic force microscopy (AFM) cantilever in single-cell force spectroscopy (SCFS), two micropipettes in dual pipette aspiration (DPA) and a microstructure fabricated using two-photon polymerization technique in single-cell adhesion micro tensile tester (SCAμTT).(d) Microplate assay is a commonly used method to monitor the force generated by individual cells when exposed to variations in substrate stiffness.Other techniques, such as micropipette aspiration, AFM and optical and magnetic tweezers, involve applying controlled displacement as the mechanical stimulus and monitoring the forces within cells or investigating changes in displacement in response to controlled forces.cells.For example, micropipette aspiration techniques have been used to investigate the microrheology of the cell nuclei [52][53][54].The monitored stress/strain response of the cellular materials will be the outcome of both active force-generating and passive mechanisms operating at those scales.From a modelling perspective, capturing active behaviour requires a mechanism of introducing change to the otherwise fixed viscoelastic models.Researchers over the years have incorporated active adaptation mechanisms by introducing active empirical mathematical models to link mechanical characteristics to the underlying biological processes without simulating the details of the underlying chemical signals.Rheological models for describing the mechanical behaviour of cellularized materials across various length scales have previously been reviewed [4,15,17,55].In this article, our primary goal is to focus on different forms of integrating activity into rheological models, emphasizing variations in their definitions.We here summarized these reported mechanisms into two broad categories.The first group integrates a force/stress (force divided by cross-sectional area)-generating element to classic viscoelastic models.This active element can be constant or time-dependent.The second group introduced mathematical methods that adjust the resting length of an elastic spring or the reference stress-free shape in two dimensions/three dimensions in classic viscoelastic models.Mathematical models in both categories are developed in conjunction with specific techniques that probe the active cell response to various types of mechanical stimuli.In some cases, the active models can be mathematically equivalent, but different parameters may offer a better link with the underlying biology.In this article, we present a selection of such models, introducing first the experimental findings and then focusing on the formulation and application of each type of active viscoelastic model.We will begin with a short introduction to the biology of active cell behaviour.

Biology of cellular active behaviour
Cell cytoskeleton mainly comprises filamentous proteins that preserve the cell structure, arrange organelles, and resist, transmit and generate forces [56,57].These proteins can be categorized into three groups: microtubules, intermediate filaments and actin filaments [58].Active force generation is realized by actin filaments.They are constructed by assembling monomeric actin.Actin-binding proteins bind to actin filaments and form different structures, such as the lamellipodium network, contractile bundles of stress fibres and the contractile network of the cell cortex [57,59].Actin filaments are engaged in active processes such as actin treadmilling and force generation by consuming the energy provided by adenosine triphosphate (ATP) hydrolysis [60][61][62].ATP molecules attach to ATP binding sites on actin monomers, and the ATP-bound actin monomers will be assembled at the plus/barbed end of the actin filaments leading to filament growth [60,63].ATP molecules will slowly hydrolyse to adenosine diphosphate (ADP) and the ADP-bound actin monomers will start to disassemble from the minus/pointed end of the actin filaments [60,64].The process of assembling and disassembling actin monomers is referred to as actin treadmilling [65].Actin filaments use the energy from ATP hydrolysis to generate protrusion forces to help cells in spreading and migration [66,67].Actin treadmilling is also crucial in endocytosis, exocytosis and phagocytosis to engulf large particles [57].
In addition, ATP hydrolysis provides the energy for contractile stress generation in the actomyosin network.Myosin motors convert chemical energy from ATP hydrolysis to mechanical energy and slide actin filaments past one another to produce force [68][69][70].This process is similar to the shortening of sarcomeres in muscle cells [37,38,64,70].The contractile forces generated from this process will then be transmitted to neighbouring cells and the ECM via cell-cell junctions and focal adhesions, respectively [66,71,72].At the cell and tissue level, contractile force production controls cellular activities, such as cell migration [67,[73][74][75], proliferation [76,77], stem cell lineage determination [78,79], tissue regeneration [80] and morphogenesis [81,82].Important for cellular mechanical characterization, contractile forces also regulate cellular response to substrate stiffness and mechanosensing [83].The active ATP-dependent processes that are at play in subcellular scales can also regulate cellular response under various loading conditions at cell and tissue scales [5,59,84,85].

Active contractile elements
Rheological models are useful tools to quantitatively analyse the results of mechanical tests.However, associating the molecular scale origin with the model parameters is challenging, in particular, where cell behaviour deviates from the response of traditional passive materials [17].Thus, adding elements that represent the active behaviour of cells can equip classic models with the ability to empirically link the mechanical response to the underlying biological processes.
To represent active cellular dynamics from intracellular contractility, contractile/force-generating elements (CEs) can be added to standard viscoelastic models.In this section, we will discuss active CEs that can mimic cell behaviour when they are exposed to different types of mechanical stimuli.We begin with the first active model, which was presented by Hill to elucidate muscle contraction and its application in studying mechanosensing in single cells.In the following section, we will explore an active fluid model that considers the transient nature of the actin network to describe characteristics of the single-cell response to the substrate stiffness.In this model framework, tension generation via myosin activity will be denoted by σ a , and regulating mechanisms of the dynamics of tension build-up will be discussed.Additionally, we will also provide some examples for both models.

Muscle contraction defined by the Hill model
Muscles function as intricate biochemical mechanisms that transform chemical energy into mechanical energy through actomyosin interactions to provide movement in our bodies.Hill proposed the first theory to describe muscle contraction back in 1938 [29].He conducted experiments using the sartorius muscle of a frog and explained the macro properties of muscles.In these experiments, Hill quantified both force generation and velocity of length changes in muscles under various loading conditions.The two endpoints of the curve representing the force-velocity relationship in muscles (figure 2a(i)) were generated in two extreme cases: zero force (isotonic condition) and zero velocity (isometric condition).Specifically, when force is kept at 0, the muscle can reach its maximum shortening velocity, V max .Similarly, when the muscle is restricted from changing its length, the maximum level of force, that is, the stall force F max , can be generated at the steady-state condition.
Hill proposed an empirical function to describe the correlation between the active force generation, F, and the shortening velocity, V, as follows [29]: where a, b and F max + a ⋅ b = c are constants specific to each muscle type.The Hill model was first discovered on frog skeletal muscles.However, later studies have shown that the dimensionless form of this model with a shape factor r = b/V max = a/F max ≈ 0.25 can be used to describe the behaviour of other muscle types, [88].Hill also proposed a phenomenological model to describe muscle mechanics.The original model consists of a contractile element (CE) that generates force and an elastic element (SE), as shown in figure 2a.CE is governed by equation (3.2) and SE affects the length and rate of change in the length of CE during contractions [29,86].Other forms of this model have also been introduced by integrating more elastic and viscous elements to account for the viscoelastic properties of the muscle and its interactions with the connective tissue surrounding muscle fibres [89].The Hill model has enabled researchers to explore the mechanics of muscles using only a few rheological parameters.However, the Hill model falls short of elucidating the underlying biological mechanisms of force generation in muscles [89][90][91].Another shortcoming of this model is that it fails to consider variations in the contractile characteristics of various fibre types within muscles and the dependence of muscle tension on the movement history [89,91].Modifications aimed at enhancing the accuracy of the Hill-type model predictions are reviewed in [91].This model is extremely versatile and not only can it be used to describe the behaviour of various types of muscles, but also it can be adapted and employed to describe force generation in non-muscle cells and their response to the physical and mechanical properties of their microenvironment [92][93][94].Mitrossilis et al. [37] studied the response of a single C2.7 myoblast cell and a 3T3 fibroblast mounted between two microplates and showed a time-dependent force generation in single cells in response to the microplate stiffness (figure 2b(iii)).The force-velocity curves for single cells have a similar shape as the force-velocity curves reported for muscle fibres.Moreover, they showed that the reduced royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231074 Hill equation (3.2), with the same shape factor r = 0.25 as for muscles, can describe the normalized shortening speed, v = V /V max versus the normalized force, f = F /F max , as shown in figure 2a(ii).In these experiments, the maximum force and velocity are measured at infinite microplate stiffness and very low stiffnesses, respectively.Therefore, showing that the coupling between force generation in cells in contact with their surroundings still follows Hill's model is quite remarkable, as it implies that structured muscle actomyosin and cytoskeletal actomyosin exhibit similar behaviours independent of the network architecture.
In addition, since the same function could explain the experimental results of different muscle types, Huxley suggested that the force-velocity relationship is generic in muscles and proposed a molecular explanation [95].This molecular model incorporates the dynamics of the interaction between actin filaments and myosin motors, that is, the number of myosin heads connected to actin filaments and the formation of temporary connections between actin and myosin heads (i.e.cross bridges).Understanding these dynamics brings valuable insight into the biological processes regulating the predictions of the Hill model.For instance, cross-bridge-type models have shown that the maximum speed of muscle shortening happens owing to the rate of myosin attachment and detachment [89,95] rather than the extent of filament overlap [96,97].This model has its limitations as well.For example, this model does not include the effect of power stroke (i.e. a crucial step in the force generation cycle in muscles where myosin heads pull the actin filaments and generate force) [89,91].
Further attempts have been made to incorporate various aspects of molecular biology with the Hill model.An example involves a model that integrated the dynamics of actin filament slippage during the process of force generation, allowing the prediction of the relationship between force generation and stiffness of micropillars [98].This model was further refined to account for the temporal evolution of the force over time by integrating an internal variable, specifically representing the progression of myosin motor stalling over time [99].Despite these improvements, this model falls short in explaining the limiting factors for the maximum force generation and maximum speed of shortening of the cells.A model that addresses these shortcomings is discussed in the next section.

Active element that includes the dynamic turnover of the actomyosin network
Étienne et al. developed an active fluid model to link the fundamental features of cellular molecular mechanics with a passive phenomenological model [38].These features encompass the transient nature of the actomyosin network, force generation via myosin motor activity and actin polymerization.The transiently cross-linked actomyosin network is modelled as a Maxwell fluid in series with an active element σ a , mimicking force generation via myosin motors (figure 2b).The constitutive equation of the network is shown as follows: Here, τ α is the characteristic time scale of cross-linker unbinding (i.e.elastic-like in short time scales (t < τ α ) and viscous over longer periods of time (t > τ α )), and E is the elastic modulus of the cell.σ a represents the maximum value of contractility or stall force that could be generated in cells, determined by two factors.First, it is influenced by the rate at which myosin motors can generate stress (1/τ myo ) and contract the cell with an elastic modulus of E. Second, it is affected by the rate of cross-linker unbinding (1/τ α ), which counteracts the increase of stress in the system.Consequently, σ a is proportional to (τ α /τ myo )E.This model effectively captures the evolution of force generation in cells leading to the establishment of tension in the steady state.Using a single dashpot (figure 2b(i)), instead of the proposed two (figure 2b(ii)), results in a similar constitutive equation.However, using two dashpots highlights the loss of force generation in the steady-state condition when the net displacement of the microplates is 0, indicating internal creep.
This model can predict the critical stiffness over which the plateau force (F p ) remains constant, the response of cells to step changes in substrate stiffness, and the rate at which force is generated across various substrate stiffnesses.Interestingly, the constitutive equation (3.3) can also be written in a similar form as the Hill model of muscle contraction, where a = E, b = 2v t + v α , and c = σ a + E v α − (LF ˙)/(2S) .Here, v = L ˙ is the speed of shortening, v t represents the rate of actin polymerization and v α = L/2τ α is the internal creep.This equation allows the examination of non-muscle cell behaviour in two extreme cases: F = 0 (stiffness of the substrate equal to 0) and v = 0 (stiffness of the substrate equal to infinity).These extreme scenarios illustrate the role of the molecular mechanisms that both govern and limit cell responses.Basically, cells initiate force generation upon attachment to the substrate.If the resistance of the external environment is lower than the force generated via myosin motors (e.g. when k = 0), cells will start to contract the microplate, which would, in turn, increase its resistance against the cell.This increased external resistance leads to a reduction in the rate of retrograde flow.Moreover, cross-linker unbinding and actin polymerization are two internal mechanisms that antagonize the rate of retrograde flow.This explanation can also be shown according to equation (3.4); when F is 0, the maximum shortening length is These two internal factors also determine the maximum force generated in the system in the case of very high stiffnesses.As demonstrated by , the maximum force in the royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231074 cell is not equal to σ a .In this scenario, actin polymerization requires extra work by myosin motors, which will be lost as a boundary creep.Furthermore, cross-linker unbinding will also result in force dissipation, leading to internal creep.An analogous active Maxwell fluid model was employed to understand the mechanical characteristics underlying the observed cellular response to local RhoA activation [87].In their experiments, Oakes et al. used optogenetic probes to recruit a cytosolic photo-recruitable protein RhoA-specific guanine exchange factor (prGEF) to the plasma membrane and activated RhoA over periods of 15 minutes.During local activation of RhoA, the fluorescent intensity of both actin and myosin II increased exponentially, plateaued and then decreased during the relaxation period (figure 2c(i)).A similar response has been observed while measuring traction forces and strain energy in cells (figure 2c(iii)).Besides, RhoA activation resulted in a sudden enhancement of traction force generation in the cell borders, whereas in the activation region traction forces did not change.
In the model used to describe this behaviour, contractile stress is introduced as an internal boundary condition, that is, l 0 < L, σ x = ± l 0 , t = ∓ σ 0 t , to represent RhoA activation in the cell area.The profile of σ 0 follows a similar pattern as the actin and myosin accumulation in the activation area (figure 2c(ii)).In addition, a two-dimensional model of stress fibres embedded in a passive viscoelastic environment was used to estimate the direction and magnitude of the actomyosin flow towards the activation region.Finally, through the application of this model, they showcased the role of Zyxin in regulating the time scale of the initial elastic behaviour.
Active CEs are employed in chemomechanical models as well.These models integrate the effect of the interplay between the mechanical characteristics of cells, external and internal mechanical stimuli and the cascade of biochemical signals to simulate cell behaviour in various scenarios.For example, a chemomechanical model illustrated the growth dynamics of cell-ECM adhesion structures and highlighted the regulatory influence of the stiffness of the nucleus and ECM [100].Another model demonstrated how the interaction between cells and ECM affects both gene expression and nuclear architecture [101].This group of models is reviewed in [102,103].

Integrating activity in viscoelastic solid models
The previous section demonstrated how active tension originating from actomyosin dynamics can be integrated within constitutive equations that relate force and rate of contraction.Étienne et al.'s model provides, in particular, a detailed description of the transient regimes leading to the establishment of steady active stress [38].These descriptions treat active materials as fluid without a reference to an intrinsic shape.It is, however, common for tissues to exhibit solid-like characteristics, with a well-defined reference shape, and the role of active stresses may then be interpreted as an apparent tension or a change in reference shape.In this section, we review these different approaches, highlighting their similarities and differences.

Constant active element
A simple way to introduce active contraction across a material is to include a stress-generating unit of constant value.Such an approach was used by Wyatt et al., when they investigated the short-time-scale response of Madin-Darby canine kidney (MDCK) monolayers to in-plane compressions [18], a process observed during morphogenetic processes [104,105] and the normal physiological function of many epithelial tissues [106,107].The MDCK monolayers were placed between two rods and compressed at different rates.Quickly compressing the monolayer to strains below a threshold level (ε ~ 33%) resulted in transient folds that disappeared in time scales of the order of seconds, while the folds created owing to higher strains were permanent.The same buckling threshold of ε ~ 33% was observed when the monolayers were compressed both rapidly and at a low rate.In addition, using actomyosin inhibitors, they have demonstrated that actomyosin activity regulates the rate of tissue flattening, the buckling threshold, pretension and the long-time-scale stiffness of the monolayers.
A simple active rheological model could reproduce the results of their studies conducted under different loading conditions (figure 3a).The model consists of a constant active element σ a , which brings the system to a tensile state even at 0 external load, in parallel to a standard linear solid model.The MDCK monolayer buckles when the stress in the monolayer approaches the compression range.Therefore, under compressive strains that would normally cause compressive stress, the model assumes that stress in the monolayer remains 0. The constitutive equation under tensile and compressive stresses is defined as , where σ a is the pretension in the monolayer, ε d is the device strain, E characterizes the long-term stiffness of the tissue, Y and η describe the short-time-scale response, ε 1 is the strain in the spring with stiffness Y satisfying ε ˙1 + ε 1 τ = ε ˙, and τ = η/Y is the characteristic time scale.Based on the model predictions for the steady state, when the applied strain is larger than the buckling threshold ε b , ε > ε b , Active rheological model that includes a constant CE, σ a , in parallel to a standard solid model to account for the prestress in cell monolayers [18].Temporal evolution of strain (i) and tissue stress (ii) as a function of device strain presented for epithelial monolayers undergoing compression at a low strain rate (0.5% s −1 ).(b) The active element proposed by Muñoz et al., which adapts its resting length in response to deformations [108].The evolution of stress over time in stress relaxation tests predicted by the Maxwell model and active models with two definitions of elastic strain is presented in a table.Stress in the models is defined using the elastic strain, that is, σ t = Eε e .Changing the definition of elastic strain can result in an equivalent evolution of stress in the active element and the Maxwell model for γ ≈ τ −1 .When the strain is very small, the difference between the evolution of stresses in all three models will be negligible.However, at larger strains, the stress in the system would be dependent on the way that the strain is defined.(c) The modified version of the active element proposed by Muñoz et al., which also considers the effect of cellular contractility [109].(i) The average curves representing the evolution of stress over time for MDCK monolayers stretched at 75% s −1 to various amplitudes of strain plotted on a semi-log scale.(ii) Demonstration of the correlation between the τ model calculated from the active rheological model and the characteristic time τ calculated from fitting the stress-time curves with an empirical function, that is, σ = Ae −t/τ t −α + B, for different loading conditions and actomyosin treatments.(d) An example of a hybrid vertex model employed to study wound healing.The top view of the epithelial tissue with ablated cells is shown along with the rheological models of nodal segments (cell-centre connections) in black and vertex segments (cell boundaries) in red [110].
royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231074 tissue strain will be equal to the device strain, ε = ε d and the stress in the tissue will follow σ = Eε d + σ a .However, when the monolayer is compressed past the model buckling threshold, stress levels fall to 0 and tissue strain plateaus at ε b = −σ a /E.These model predictions are consistent with the experimental results shown in figure 3a(i) and (ii).The model can also predict the response of the monolayer to a step of compressive strain of different magnitudes.Besides, the model provides a simple way to capture how treatments affecting actomyosin control the amount of active tension revealing that the buckling threshold observed in the experiments agrees with model predictions −σ a /E in all conditions.

Active behaviour modelled as change in the resting length of a spring
The previous approach focused on how active stress, combined with external perturbations, would set the mechanical evolution of a viscoelastic tissue.The Maxwell branch and its dashpot, in particular, account for the remodelling and plasticity of the material.Another way to represent this is through the explicit evolution of the rest shape, or stress-free state, of the material.Both remodelling and tensioning can be accounted for through the evolution of the resting length of cells and tissues over time, as proposed by Muñoz et al. [108].We first describe how a Maxwell-like behaviour can emerge from this strategy and then how activity can be added to this class of models.

Model proposed by Muñoz et al.
Muñoz et al. proposed a model to account for the influence of cytoskeleton activity on cell shape by representing plasticity as a change in the rest length of a spring [108].A material's constitutive equation therefore takes the form of a relationship between the rate of change of the rest length and the stress (or equivalently some metric of elastic strain) in the material.In this model, the current resting length of the material (filament, cell or tissue) L, that is, the total length of the material when no external load is applied to it, is proportional to the elastic strain ε e .The rate of changing the resting length under strain is defined as where γ is the remodelling rate of the network defined as the network resistance to adjusting its configuration to the applied deformation, ε e is the current elastic strain, ε e = l − L L and l is the current total length of the network.According to Muñoz's model, the definition of current elastic strain is different from the apparent strain ε = l − L 0 L 0 , where L 0 is the initial length and the resting length of the network.For small deformations, Muñoz's approach is mathematically equivalent to a linear Maxwell model with a characteristic time of τ = η E [108].However, at large deformation, the models differ, leading to rather complex relaxation dynamics for Muñoz's model (see figure 3b).This difference results from the particular definitions of strains and the resulting nonlinearities emerging from them.For instance, using a logarithmic form for the elastic strain ε e = ln l L (i.e. using the true strain definition rather than the engineering strain), Muñoz's model would match the Maxwell model up to large deformations, with a relaxation time scale independent of the strain amplitude, as shown in the table in figure 3b.At a mathematical level, controlling the rest length of a spring or having a dashpot in series is therefore largely equivalent.However, introducing a dynamic rest length enables a slightly different interpretation of the physiological mechanisms [108].The dashpot element in the Maxwell model is often associated with remodelling the system but could also account for the resistance of the cytoplasmic fluid η to the applied strain rate and dissipate power.On the other hand, for Muñoz's active element, the inelastic part of the external power will be used to overcome the resistance of the cytoskeleton filaments γ to adapt to the new configuration imposed by the external strain.The active model proposed by Muñoz has been generalized to two-dimensional/three-dimensional continuum models [111] and also integrated into discrete models such as cell-centred [112], vertex [113] and cell-centred/vertex hybrid [109,110] approaches.By incorporating a porosity parameter representing the density of polymers in the cell cytoskeleton, the continuum model proposed by Asadipour et al. can also replicate the immediate fluidization in cells in response to transient strains and the subsequent gradual stiffening 10 royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231074 [111].These adaptations have facilitated the study of epithelial tissue behaviour in both two and three dimensions.A modelling approach similar to the Muñoz model proposed by Esfahani et al. could also demonstrate stiffening in response to high strain rates applied to epithelial cell pairs [49].A few of the modifications made to the Muñoz model are presented in the following sections.

Active element considering cellular contractility
Mosaffa et al. [109] modified the evolution law of the resting length of the material proposed by Muñoz et al. [108] by introducing a contractility parameter ε c to account for the inherent contractility of the cells, In this model, when the elastic strain ε e reaches to ε c , the resting length will not change any more, and as previously stated [108], when ε c is zero, the model behaves similarly to the Maxwell model.Mosaffa et al. implemented the modified active element in a hybrid cell-centred/vertex model where cells interact through both cell centrs, presented by nodes and cell-cell junctions, presented by the connection between vertices [109].This model could successfully simulate tissue extension and wound healing.
Khalilgharibi et al. have used a similar approach for fitting the results of stress relaxation tests conducted on MDCK monolayers [5].Their studies have shown that stress in the MDCK monolayer increases promptly after strain application.Then, the stress will gradually relax along with an increase in the monolayer length, which is regulated by actomyosin activity.Moreover, they have noticed a strain-dependent characteristic time τ for monolayers stretched at 75% s −1 strain rate (figure 3c(i)), which cannot be explained by standard linear viscoelastic models.Therefore, they proposed a model that consists of an elastic spring in parallel to an active element that sustains a constant pre-strain ε c and changes its resting length L t to relax the imposed stress and return its strain to ε c .The changing of the resting length of the monolayer L t in response to an applied strain ε 0 is defined as Here, ε e is the effective strain ε e t = l m − L t /L t with l m representing the actual length of the monolayer after applying the deformation and γ is the rate of changing the resting length.The characteristic time predicted by this model, τ model = ε 0 / γ 1 + ε 0 , increases with the applied strain, which is consistent with their experimental observations, as shown in (figure 3c(ii)).Another example of the active element that accounts for cellular contractility involves integrating this model into a hybrid two-dimensional cell-centred/vertex model to analyse the wound healing process [110].In this model, the vertex segments (cell boundaries) and nodal segments (cell-centre connections) connect the apical and basal sides.Nodal segments are characterized using an elastic spring, and the behaviour of vertex segments is described by an active model consisting of two branches in parallel (figure 3d).The first branch is an active element that accounts for the changes in the resting length of the vertex following equation (4.3).The second branch is an elastic spring with an additional time-varying contractility parameter Υ ^c, which accounts for the effect of high actin concentration at the wound edge (i.e.purse string contractility) and increases the stress in the elastic spring.The outcomes of their simulations have demonstrated the regulatory mechanism of both purse string contractility and tissue contractility on wound healing speed.

Active element considering a time delay in the active rest length changes
One of the potential reasons for the oscillatory response observed in tissues during various processes, such as morphogenesis, could be the delay between the signal and the response, as stated by Muñoz et al. [113].These delays in the responses can be owing to the distance between the sender and receiver of the biochemical signals or the time necessary for signal processing [113,114].To study this phenomenon, they modified the active element that was previously proposed by Muñoz  τ oscil is the time beyond which the rest length of the element oscillates, e is the exponential constant, and for time scales above τ stab the value of the rest length is unstable, and its oscillation amplitude will rise over time.These values might trigger oscillations during embryogenesis.
In addition, delays can also be dependent on the apparent size of the element l t .Muñoz et al. implemented the effect of the size-dependent delays, τ t = λl t , into a vertex model to analyse the oscillation in the cellular area in biological tissues [113].In the vertex model, the rest length of the nodal elements was maintained constant, and changes in the resting length of the vertex elements were defined using equation (4.5).For constant delays, oscillations in the cell area were periodic and synchronous.However, for size-dependent delays, oscillations started to get increasingly out of sync.It is argued that this model demonstrated the role of delay in the mechanical response in inducing oscillations even in the absence of external sources.

Summary and future perspectives
The active viscoelastic models outlined here are capable of capturing numerous aspects of cell behaviour at multiple spatial and temporal scales with a small number of model parameters without considering the details of the structural components and dynamics of cell-cell and cell-ECM adhesion sites.These active models provided insight into the results of experiments and predicted the system behaviour in other arbitrary conditions.Additionally, despite the difficulty in establishing a clear connection between biological processes and model parameters, researchers have used drug treatments or targeted mutations to demonstrate correlations between model parameters like Young's modulus, viscosity or active pretension and biological processes such as actomyosin activity, even molecules that regulate these processes.
Each active viscoelastic model is described by a constitutive equation that represents a particular cellular behaviour, such as sensing changes in the substrate stiffness, and response to strains at different magnitudes and rates in cell doublets and cell monolayers.Therefore, the existence and use of a generalized model that can be employed to describe and predict the response of cells in different scenarios is still an open question.Integrating and bridging the gap between phenomenological and biophysical models is an important step to improve our understanding of these systems.For example, a molecular model of the actomyosin cortex inspired a phenomenological model for cell-scale mechanosensing [38], and a phenomenological model that included the role of actin polymerization in changing the resting length of the material was incorporated in vertex models to study different aspects of tissue dynamics [113].Consequently, a comprehensive phenomenological model that captures all the significant facets of rheological data might also enrich the findings of biophysical models and allow us to improve the precision of simulations of cell activity.
Numerical models are an excellent research tool to complement, analyse and interpret experimental data in the field of cell mechanics.Numerous models with varying degrees of complexity and details of the structural elements involved in the observed phenomena have been presented over the years.Power law [15], fractional (reviewed in [17]), viscoelastic and active viscoelastic models are included in this category, where the effects of subcellular microstructures on cell rheology are represented by model parameters.For example, virtual cell (VCell) is a powerful model that includes details of the nucleus, cytoskeleton, cytoplasm and chromatin fibres [115].This level of detail is computationally expensive and might not always be necessary.In other words, based on the research question, the phenomenon of interest, length and time scales and characteristics of the relevant microstructures can be incorporated into biophysical models.For example, at the molecular scale, chemomechanical [103,116] and molecular clutch models [117][118][119] are used to study cell-cell and cell-ECM adhesions, at the cell-scale statistical approaches can be employed to study cell mechanics [120][121][122], and cellular Potts (CPM) [123][124][125], vertex [126,127] and self-propelled Voronoi [128][129][130] models are introduced to study the mechanical behaviour of epithelial monolayers in two dimensions and three dimensions.The level of detail can be further reduced by describing the outcomes of experiments using phenomenological models before examining the underlying mechanisms.

Figure 1 .
Figure 1.Experimental techniques used to probe cell mechanics in various spatial scales.(a) Mechanical testing machine used to study the response of tissue samples in compression and tension tests.(b) The recent technique developed for investigating the mechanical response of cell layers to uniaxial stretch and compression.(c)In cell pair studies, mechanical stimulus is applied through various methods such as a single cell adhered to a tipless atomic force microscopy (AFM) cantilever in single-cell force spectroscopy (SCFS), two micropipettes in dual pipette aspiration (DPA) and a microstructure fabricated using two-photon polymerization technique in single-cell adhesion micro tensile tester (SCAμTT).(d) Microplate assay is a commonly used method to monitor the force generated by individual cells when exposed to variations in substrate stiffness.Other techniques, such as micropipette aspiration, AFM and optical and magnetic tweezers, involve applying controlled displacement as the mechanical stimulus and monitoring the forces within cells or investigating changes in displacement in response to controlled forces.

Figure 2 .
Figure 2. Models of active contractility relating force/tension generation and the rate of contraction.(a) Model proposed by Hill[29] to describe muscle contraction[86].CE serves as a force-generating contractile element and SE is an elastic spring that affects the dynamics of contraction.(i) Force-velocity curve obtained in Hill's experiments.(ii) Normalized speed of shortening v = V /V max and normalized force generation f = F /F max in single cells in response to various substrate stiffnesses can be represented by the reduced Hill equation (3.2)[37].(b) Active fluid model proposed to study mechanotransduction in cells[38].(i) This model accounts for the transiently cross-linked nature of the actomyosin network and force generation via myosin motor activity.(ii) Equivalent of the model presented in (i), which has the same constitutive equation.The net effect of two dashpots is similar to one dashpot.(iii) Evolution of traction force generation over time for low (grey) and high (black) substrate stiffness.The two main characteristics of force generation-time curves, that is, the plateau force F p and the rate of force generation dF /dt are dependent on the substrate stiffness[37].(iv) The plateau force has a linear relationship with substrate stiffnesses when the stiffness is below 60 nN μm −1 and for higher values of stiffness, F p saturates[37].(c) Active Maxwell model that assumes a time-dependent increase in cellular contractile force with the same profile as actin and myosin accumulation in the RhoA-activated areas[87].(i) Changes in the local intensity of actin and myosin in the regions of RhoA activation.(ii) The plateauing exponential function used in the model to represent the stress profile during activation of RhoA.(iii) Experimental (black) and theoretical (red) evolution of strain energy during intervals of RhoA activation and relaxation.

Figure 3 .
Figure 3. Active models for solid-like materials.(a).Active rheological model that includes a constant CE, σ a , in parallel to a standard solid model to account for the prestress in cell monolayers[18].Temporal evolution of strain (i) and tissue stress (ii) as a function of device strain presented for epithelial monolayers undergoing compression at a low strain rate (0.5% s −1 ).(b) The active element proposed by Muñoz et al., which adapts its resting length in response to deformations[108].The evolution of stress over time in stress relaxation tests predicted by the Maxwell model and active models with two definitions of elastic strain is presented in a table.Stress in the models is defined using the elastic strain, that is, σ t = Eε e .Changing the definition of elastic strain can result in an equivalent evolution of stress in the active element and the Maxwell model for γ ≈ τ −1 .When the strain is very small, the difference between the evolution of stresses in all three models will be negligible.However, at larger strains, the stress in the system would be dependent on the way that the strain is defined.(c) The modified version of the active element proposed by Muñoz et al., which also considers the effect of cellular contractility[109].(i) The average curves representing the evolution of stress over time for MDCK monolayers stretched at 75% s −1 to various amplitudes of strain plotted on a semi-log scale.(ii) Demonstration of the correlation between the et al. by 11 royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.11: 231074 considering the effect of a time delay between the mechanical signals and the active rest length changes [108], (4.5) L ˙= γ l t − τ − L t − τ .Analysing the stability of the delay differential equation resulted in the limits of oscillation and stability as follows: